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ALGORITHM  FOR  THE  INVERSE  OF 
A  HERMITIAN  TOEPILTZ  MATRIX 


INTRODUCTION 

The  efficient  inversion  of  a  given  matrix  and  the  related  problem  of  solving  a  system  of  linear 
equations  has  been  a  subject  of  intense  study  for  many  years.  The  literature  on  this  subject  is  so  vast 
that  no  survey  can  be  exhaustive.  For  example,  a  tentative  classification  and  bibliography  on  solving 
systems  of  linear  equations  written  by  Forsythe  ll]  contains  over  400  titles.  An  excellent  handbook  on 
the  various  numerical  methods  of  matrix  inversion  and  the  solution  of  linear  equations  has  been  writ¬ 
ten  by  Westlake  (2).  Different  methods  are  compared  based  on  such  measures  of  effectiveness  as 
speed,  storage  requirements,  and  convergence  rales  if  applicable. 

Numerical  methods  for  matrix  inversion  and  the  related  problem  of  solving  a  system  of  linear 
equations  can  be  divided  into  two  classes:  the  direct  methods  and  the  indirect  (iterative)  methods. 
Direct  methods  such  as  Cramer’s  rule  13],  Gaussian  elimination  [3],  and  orthogonalization  [3-4]  yield 
an  exact  solution  after  a  finite  number  of  operations  if  there  is  no  roundoff  error.  Iterative  methods  on 
the  other  hand  such  as  gradient  methods  [4],  the  back  and  forth  Seidel  [4],  and  successive  overrelation 
[5],  begin  with  an  approximate  solution  and  obtain  an  improved  solution  with  each  step  of  the  iteration. 
The  accuracy  of  the  solution  depends  on  the  number  of  iterations  performed. 

For  most  direct  methods  of  matrix  inversion,  the  number  of  arithmetic  operations  is  proportional 
to  M 3  where  M  is  the  row  or  column  dimension  of  the  given  square  matrix.  For  iterative  methods,  the 
number  of  operations  per  iteration  is  proportional  to  M1.  In  general,  the  speed  of  an  algorithm  if  there 
is  no  parallel  processing  is  proportional  to  the  number  of  arithmetic  operations  so  that  this  measure  can 
be  used  to  evaluate  the  performance  of  a  given  algorithm. 

A  direct  procedure  for  finding  the  solution  of  simultaneous  linear  equations  where  the  multiplying 
matrix  is  Toepiltz  was  developed  by  Levinson  and  presented  in  Norbert  Weiner’s  book.  Extrapolation , 
Interpolation,  and  Smoothing  of  Stationary  Time  Series  16].  This  algorithm  takes  advantage  of  the  Toepiltz 
form  to  reduce  the  number  of  arithmetic  operations  to  be  proportional  to  M2.  This  algorithm  has  been 
used  by  Burg  [7]  to  estimate  line  spectra  in  a  methodology  commonly  called  maximum  entropy  spec¬ 
trum  analysis  (MESA). 

This  report  presents  a  new  direct  method  for  finding  the  inverse  of  an  M  x  M  hermitian  Toepiltz 
matrix.  An  M  x  M  hermitian  Toepiltz  matrix,  H,  has  the  form 
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where  h0  is  always  real  and  *  indicates  the  complex  conjugate.  Note  that  it  is  only  necessary  to  specify 
the  elements  of  the  first  column  of  a  hermitian  Toepiltz  matrix  in  order  to  define  the  entire  matrix. 
Therefore,  we  introduce  the  shortened  notation:  if  H  is  an  M  x  M  hermitian  Toepiltz  matrix,  then  we 
write 


H  “  ((Ao>  h\,  /ij . |))  (2) 

where  hk,  k  —  0,  1 . M  -  1  are  the  elements  of  the  first  column  of  H. 

We  will  take  advantage  of  the  form  of  a  hermitian  Toepiltz  matrix  and  develop  new  direct 
methods  for  the  solution  of  simultaneous  linear  equations  and  the  matrix  inverse.  The  basis  of  these 
related  algorithms  lies  in  discrete  Fouries  series  theory.  Efficient  algorithmic  procedures  are  presented 
which  use  the  theory  of  the  preceding  sections  to  find  the  matrix  inverse  and  the  solution  of  simultane¬ 
ous  linear  equations  respectively.  We  then  discuss  the  software  implementation  of  the  matrix  inversion 
algorithm. 

SOLUTION  OF  SIMULTANEOUS  LINEAR  EQUATIONS 

In  this  section,  we  will  develop  an  algorithm  for  solving  for  the  unknowns  of  a  system  of  M 
independent  linear  equations.  Using  this  algorithm,  we  will  see  in  the  next  section  that  an  algorithm 
for  obtaining  the  inverse  of  a  given  hermitian  Toepiltz  matrix  can  be  derived. 

Consider  the  vector  equation 

m  =  r  (3i 


where  //  is  an  M  x  M  nonsingular  hermitian  Toepiltz  matrix,  c  is  an  ,W  x  1  known  vector,  and  x  is  an 
M  x  1  unknown  vector.  We  desire  to  find  x  We  use  the  following  approach.  Let  us  define  a  system 
of  iV  =  2 M  -  1  independent  linear  equations  as 
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where  xM  is  an  M  x  1  unknown  vector,  cM  =  c,  xM-i  is  an  (M  -  1)  x  1  unknown  vector,  0  is  a  zero 
filled  (M  -  1)  x  1  vector,  and 
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We  call  a  matrix  defined  by  the  form  seen  in  Eq.  (5)  as  an  up-down  hermitian  Toepiltz  matrix 
(UDHTM)  because  the  subscripts  of  hk  seen  in  Eq.  (5)  increase  and  then  decrease.  We  also  assume 
that  PKt  is  nonsingular.  The  matrix.  Pm,  as  seen  in  Eqs.  (4)  and  (S)  is  partitioned  as  follows:  />nM>  is  a 
M  x  M  matrix,  Pt2l)  is  an  M  x  (M  -  1)  matrix,  /*2tw>  >s  an  (M  —  1)  x  M  matrix,  and  P^  is  an 
(M  -  1)  x  (M  -  1)  matrix.  In  addition,  we  can  also  show  that  />i(jM)  -  H  and  thus  is  hermitian  Toe¬ 
piltz  and  that  P^1  is  also  hermitian  Toepiltz.  In  fact 

Pi}1)  ”  (.(ho-  h\,  hi.  •••  .  hui-i)).  (6) 
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Note  that  the  system  of  equations  defined  by  Eq.  (4)  contains  (2 M  —  1)  unknowns  and  has  a  unique 
solution  if  P,( xM)  is  nonsingular.  Also  note  from  Eq.  (4)  that  if  Pttfi>  -  H  and  cM  -  F,  then  xM  -  x. 
Thus,  if  we  solve  for  the  unknowns  in  Eq.  (4),  we  have  also  solved  for  x  in  Eq.  (3). 

Let  us  rewrite  Eq.  (4)  as 


where 
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(8) 


such  that  C?nM)  is  an  M  x  M  matrix,  Qf 2W)  is  an  M  x  (M  -  1)  matrix,  OjiM>  is  an  (M  -  1)  x  M 
matrix,  and  Q 22  '  is  an  (M  -  1)  x  (.M  -  1)  matrix.  We  show  in  Appendix  A  that 

QM=jjF;\'F,  (9) 

where  FlS  is  the  Mh  order  discrete  Fourier  series  (DFS)  matrix  defined  by  Eq  (A6),  and  A  is  a  diago¬ 
nal  matrix  whose  element,  Kkk ,  consists  of  the  Mh  order  DFS  of  the  sequence  |/i0,  /» 1 ,  ...  ,  hM.  1, 
hi,. |,  hi,. 2,  ■  ■  ,  hl\  (the  Mh  order  DFS  is  defined  by  Eq.  A4).  In  fact,  if 

U.  si . **-■)  -  DFS  {*0.  h  1 . hM.\.  h'M.x . /i,*)  (10) 

where  (s0.  S| . sN.x )  is  the  sequence  that  results  by  finding  the  DFS  of  the  sequence  (/j0,  h . . 

hit- . . >0.  ‘hen 

A**  —  j*_i;  A  —  1,  2 . N.  (11) 

It  is  also  shown  in  Appendix  A  that  if  PM  is  a  UDHTM  then  Pi,'  or  Qu  is  also  a  UDHTM.  Thus  we 
see  that  Qu  can  be  written 

Qm-  <(«o.  Si.  1.  Sw-i . *r»-  D2) 

Hence,  it  is  seen  that  Q {f,)  and  Qjf*'  are  hermitian  Toepiltz  matrices  with 

air  -  «8o.  «i . *«-i»  (i3a> 

and 

GW'  -  ((*0.  Ml . *M- 2».  03b) 

Now,  let  us  rewrite  Eq.  (7)  in  the  equivalent  form  as 

-  <?,r  +  Q\¥ '  x„-i  (14a) 

0-  c*  +  QW}xk-\.  d4b) 
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Equation  (14b)  can  be  rewritten  as 


p<r.) 


xu- 1 


where  we  have  defined 

-  (?2r‘;  cu-t--Qircu. 


(14c) 

(15) 


If  we  had  a  solution  for  xw_|  in  Eq.  (14a),  we  could  find  xu.  To  find  xw_i,  we  use  Eq.  (14c).  How¬ 
ever,  is  a  (M  —  1)  x  (M  -  1)  hermitian  Toepiltz  matrix  and  cM_i  is  a  derivable  (M  -  1)  x  1 

vector.  Thus,  we  have  reduced  the  order  of  the  problem  from  finding  an  unknown  M  x  1  vector,  xM, 
to  finding  an  unknown  (M  -  1)  x  1  vector,  xM~u  whose  multiplying  matrix  is  also  hermitian  Toepiltz. 
Hence,  the  above  procedure  is  reiterative  and  must  be  repeated  M  —  1  times  with  the  assumption  that 
Pn*  k  —  2 . Mare  nonsingular.  On  the  M  -  1  iteration,  the  equations  have  the  form 


*2  -  QiV  c2  +  q,t  x, 

-  ~  Qi?*  Cj 

x' "  QW 


(16a) 

(16b) 


Note  that  X|  and  Q 2<|)  are  now  scalars.  Thus,  no  matrix  inversion  of  @22*  is  necessary  (it  is  assumed 
022*  ^O).  Therefore,  X\  is  known  and  x2  (a  2  x  1  unknown  vector)  can  be  found  by  using  Eq.  (16a). 
In  general,  the  unknowns,  xk,  k  =>  2,  3 . M  can  be  obtained  by  using  the  forward  reiterative  for¬ 

mula 

**  -  (?n  1  ck  +  (){$'  xk.,.  (17) 

The  constant  k  x  1  vector,  ck,  k  =  2,  3 . M  can  be  obtained  by  using  the  backward  reiterative  for¬ 

mula 

c*_,  -  -  QiV  ck\  k-M.M-  1 . 2  (18) 


with  the  final  condition  that  cM  =  c.  In  the  discussion  of  software  algorithm  for  matrix  inversion,  we 
discuss  how  to  obtain  the  matrix,  Qk,  k  =  M,  M  -  1,  ....  2. 

TOEPILTZ  MATRIX  INVERSION  ALGORITHM 


We  can  use  the  algorithm  for  finding  the  unknowns  of  a  system  of  linear  equations  discussed  in 
the  preceding  section  to  obtain  the  inverse  of  a  given  hermitian  Toepiltz  matrix.  Let  us  define  the 
(xk  matrix,  ft*,  such  that 


n*  A  I#>n,rl  ;  *  -  2,  3 . M 


(19) 


and 


ft,  - 


1 

oir 


(20) 


Note  from  Eqs.  (3),  (4),  and  (19)  that  ft  w  -  H~\  Now  for  k  -  M,  Eq.  (14c)  implies  that 

*m  ”  ftjw  (21) 


Equations  (21)  and  (15)  imply  that 

*m-\  ”  ft«-ic«-i  “  ~  ft  m—\  Qi\,)  cm  m  ~  ftjM-i  Qm-  1  Qir  c.  (22) 
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Therefore,  if  we  substitute  Eqs.  (21)  and  (22)  into  Eq.  (17),  we  obtain 

C  -  Qir  ft  u Qir  C  (23) 

-  (Qif*  ~  Q i'i"* 

Because  c  is  an  arbitrary  M  x  1  vector,  Eq.  (23)  implies  the  following  formula: 

ft*#-  -  Q\?'*u-\QkV'-  (24a) 

Similarly,  we  can  find  a  formula  for  ft  M  -  •  We  do  this  by  choosing  a  new  arbitrary 

vector,  c,  of  length  M- 1 .  and  initiate  solving  a  system  of  M- 1  simultaneous  equations  as  we  did  in 
the  preceding  section.  Using  equations  similar  to  Eqs.  (21)  to  (23),  we  would  derive  an  equation 
exactly  like  Eq.  (24a)  except  that  the  index  is  M-l.  Hence  it  is  possible  to  write  a  reiterative  formula 

ft*-  O/f*  -  QiV  ft*-.  Qti'  (24b) 

with  k  =  2,  3 . M,  and  with  ft,  given  by  Eq.  (20).  Thus  if  we  reiterate  Eq.  (24b)  M  -  1  times, 

we  obtain  H~x  =  ftM. 

SOFTWARE  ALGORITHM  FOR  MATRIX  INVERSION 

In  this  section,  we  present  an  efficient  procedure  for  obtaining  the  inverse  of  a  hermitian  Toepiltz 
by  using  the  methodology  described  in  the  preceding  sections.  To  begin  with,  it  is  seen  from  Eqs.  (20) 

and  (24b)  that  all  that  is  necessary  for  computing  the  11  k  matrices,  k  =  1,2 . M  are  the  Qk 

matrices.  The  partitions,  Q|*\  On',  and  (?22 '  can  be  obtained  easily  from  Qk  Now  Qk  is  a 

UDHTM,  so  that  all  that  is  necessary  to  completely  specify  it  is  the  first  column  of  the  matrix  (actually 
because  of  the  up-down  properly,  just  the  first  M  elements  of  the  first  column  are  needed).  The  matrix 
Qk  can  be  found  by  using  the  formula 

Qk  ”  —77  Fu-t  A*1  Fm_,;  k  =  2.  3 . M  (25) 

where  's  tbe  (2*  ~  D  order  DFS  matrix  defined  by  Eq.  (A6)  and  A*  is  a  diagonal  matrix.  The 
diagonal  element,  X/5*\  /  -  1,  2,  . . .  ,  2k  -  1  is  found  as  follows.  If 

Qk+i  “  ((SoU+,).  W*+l) . £*U+1>.  . £i<*+,)‘))  (26) 

and 

UoU).  sik) . still 2}  -  DFS  (#o(*+U.  . . (27) 

then 

s/-V  ;  1-1.2 . 2k  -  1.  (28) 

Now  in  order  to  generate  all  Qk,  A  w  must  be  known.  However,  the  matrix,  \M,  can  be  obtained 
by  using  the  elements  that  define  H and  Eqs.  (10)  and  (11). 

To  evaluate  Qk,  seen  in  Eq.  (25),  it  is  not  necessary  to  perform  all  of  the  matrix  operations  indi¬ 
cated  by  this  equation.  In  fact,  it  is  straight  forward  to  show  (see  Appendix  A)  that 

. iH- 

_ i _ 1  _ 1 _ 

(2k  -  !)A/f>  ’  (2k  -  DA ,<f> .  (2k  -  1)A &>_  „(M  _  „ 
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Thus,  based  upon  the  preceding  discussion,  we  present  the  following  algorithmic  procedure  for 
finding  H~u. 


A. 

Set/c-  M-  1  h„  /-0,  1,  .. 

.  ,  M  -  1. 

B. 

Calculate  {s^**,  s\k) . su-i)  by  using  Eq.  (27). 

C. 

Calculate  {g<j*’<  g{*K  -  -  -  •  ftk-i)  by  using  Eqs.  (28)  and  (29) 

D. 

Store  (go'*’,  g/** . g*-ii. 

E. 

*  -  k  -  1. 

F. 

Go  to  B  if  *  >  1 . 

G. 

Set  ft  i  -  — y  (note  go21  -  Q221)  and  k 
go 

-  2. 

H. 

Set  Qk  -  ((go**,  g,U),  . . .  ,  g**{,  g*-i*,  •  • 

..g,u,*». 

I. 

Construct  partitions:  Gi* \  0n\  QiV- 

J. 

Otf'  -  QiV  o*_,  ftj'f*. 

K. 

k  =  k  +  1. 

L. 

Go  to  H  if  k  <  M. 

M. 

The  algorithm  can  be  divided  into  two  parts:  the  first  part  (steps  A-F)  consists  of  finding  the  elements 
of  Qk ,  k  =  2,  3,  . . .  ,  M,  and  the  second  part  (steps  G-M)  calculates  through  a  reiterative  formula  (step 
J)  the  ft*  matrices. 

It  can  be  shown  that  gbk)  and  s/*\  /  =  0,  1,  . . .  ,  2k  —  2  are  always  real.  Because  of  computa¬ 
tional  errors,  however,  these  values  may  have  a  small  imaginary  part.  It  was  found  that  the  accuracy  of 
the  matrix  inverse,  H~\  improved  if  only  the  real  part  of  the  computed  g0U)  or  s/*1  was  used  in 
succeeding  steps  of  the  algorithm. 

A  Fortran  computer  program  listing  that  implements  the  matrix  inversion  algorithm  is  given  in 
Appendix  B. 

SOFTWARE  ALGORITHM  FOR  SOLUTION  OF  SIMULTANEOUS  EQUATIONS 

Similar  to  the  preceding  section,  we  present  an  algorithmic  procedure  for  finding  the  solution  of  a 
system  of  simultaneous  linear  equations  as  given  by  Eq.  (3)  as  follows: 

A.  Set  k  -  M,  g/M+,)  0,  1,  2 M  -  1,  cM  -  c. 

B.  Calculate  {so(*\  s{k\  <  *2*-ii  by  using  Eq.  (27). 

C.  Calculate  {go*’*  g|<u,  •  •  •  .  g*-!}  by  using  Eqs.  (28)  and  (29). 
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D.  Store  !gd*\  gi*’ . &-!) 

E.  Construct  partition  Qjf*  by  using  {go*',  g{*> . g£-\). 

F  i  -  -  C?#’  c*;  store  c*_,. 

G.  k  -  k  -  1. 

H.  Go  to  B  if  /c  >  1 . 

I.  Set  x,  —  C|/foJ)  (note  x1(  ?|  are  scalars)  and  k  -  2. 

J-  Set  Qk  -  «*}*>.  g}*\  ....  git\,  g£k-V . «,“>*)). 


K. 

Construct  partitions: 

Q,(}\  0,V 

L. 

xk  Vll  ck  +  (e  12 

xk- 1- 

M. 

k  -  k  +  1. 

N. 

Go  to  J  if  k  <  M. 

O. 

1 

IX 

IMPLEMENTATION  OF  THE  MATRIX  INVERSION  ALGORITHM 

The  value  of  any  algorithm  that  is  used  as  a  Cuir.f  uter  library  subroutine  is  determined  by  such 
measures  as  speed,  the  amount  of  computer  memory  needed,  and  the  amount  of  hardware  necessary  to 
implement  the  algorithm.  The  last  two  measures  can  sometimes  be  traded-oflf  to  obtain  faster  speeds 

For  the  matrix  inversion  algorithm,  the  amount  of  memory  (double  words  for  a  complex  number) 
needed  is  at  most  M2.  To  see  this,  we  observe  from  steps  A-F  that  it  is  necessary  to  store 
M  (M  —  1 )/  2  complex  numbers.  For  steps  G-M,  it  is  necessary  to  store  at  most  M2/  2  complex 
numbers.  This  results  because  it  can  be  shown  that  if  (lt  =  (fi,S,V),  *  “  1,  2,  ....  M,  then 

ft  mV  ”  ft  of-m  +  nu-B+n;  m.n  -  1.  2 . k.  (30) 

Therefore,  only  half  of  the  elements  of  the  fl  k  matrix  need  to  be  stored.  Since  k  ^  M.  this  number  is 
ai  most  M2/2.  Hence,  it  follows  that  the  maximum  memory  needed  for  steps  A-M  is  M2.  Storage 
requirements  for  most  matrix  inversion  algorithms  are  of  the  order,  M 2  12].  Thus  there  is  no  advan¬ 
tage  in  eliminating  memory  by  using  the  matrix  inversion  algorithm  presented  in  this  report. 

A  good  indication  of  the  speed  of  an  algorithm  is  the  number  of  multiplications  (As)  that  are 
necessary  to  perform  the  algorithm.  Multiplications  and  divisions  that  are  implemented  digitally  are 
generally  much  slower  operations  than  the  addition,  subtraction,  loading,  and  storing  operations  and 
hence  may  account  for  the  greater  portion  of  the  processing  time.  For  steps  A-F  of  the  matrix  inver¬ 
sion  algorithm,  the  approximate  number  of  As  is  2A/V3  and  for  steps  G-M  the  approximate  number  of 
As  is  M4/ 4.  Hence,  the  total  number  of  As  is  of  the  order  of  M4.  For  most  direct  methods  of  matrix 
inversion  the  number  of  As  is  of  the  order  of  M3  12],  Thus  it  is  seen  that  the  algorithm  presented  in 
this  report  is  comparatively  slow  at  least  when  implemented  in  pure  software. 

There  are  two  other  disadvantages  associated  with  this  matrix  inversion  algorithm.  First,  if  the 
given  hermitian  Toepiltz  matrix  H  is  singular,  the  algorithm  does  not  indicate  this.  Second,  if  H  is 
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nonsingular,  the  intermediate  UDHTMs  employed  in  the  algorithm  may  be  singular.  In  this  ease,  the 
algorithm  fails.  It  is  possible  to  determine  if  an  intermediate  UDHTM  is  singular  by  noting  whether 
any  of  the  values  of  klh  /  -  1,  2,  . . .  ,  2k  -  1,  calculated  in  Eq.  (28)  are  zero.  If  any  of  these  values 
are  zero,  then  the  given  UDHTM  is  singular  and  the  algorithm  fails. 

SUMMARY  AND  CONCLUSIONS 

A  new  method  for  obtaining  the  inverse  of  a  hermitian  Toepiltz  matrix  was  presented.  In  addi¬ 
tion,  a  related  technique  for  finding  the  solution  of  the  system  of  linear  equations,  Hx  =  c,  where  H  is 
a  hermitian  Toepiltz  matrix,  was  developed.  Efficient  algorithmic  procedures  for  both  of  these  methods 
were  listed. 
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Appendix  A 

INVERSE  OF  AN  DP-DOWN  HERMITIAN  TOEP1LTZ  MATRIX 

In  this  appendix,  we  derive  the  inverse  of  a  nonsingular  up-down  hermitian  Toepiltz  matrix 
(UDHTM).  Let  A  be  a  N  x  \  UDHTM  such  that 

A  =  ((<?()•  ai‘  a2 . aM-\ ■  aM- 1.  aM-2 . a\  ))  (Al) 

where  a0  is  real  and  N  =  2M  -  1. 

The  methodology  of  finding  A~ 1  is  embedded  in  discrete-Fourier-series  (DFS)  analysis.  The  DFS 

periodic  convolution  theorem  [Al]  states  that  if  xik),  v(k )  and  z < A > ,  k  =  ...  —2,  —1,  0,  1,  2 . 

are  periodic  sequences  with  a  period  equal  to  N  and 

2(ti)=  ’['rW.vlii-m),  ( A2) 

m—0 

then 

Z(k)=X(k)Uk)  (A3) 

where  Xik),  Yik),  and  Z(k)  are  the  /Vth  order  DFSs  of  xin).  yin),  and  :in)  respectively.  Recall 
that  a  DFS  is  defined  by  the  mapping  of  a  sequence,  u  in),  of  length  N  into  a  sequence.  Uik),  through 
the  transformation 

Uik)  =  X  **",%"•  k-  0.  1.  2 . N  -  1  <A4) 

n=0 

where  W\  =  exp  \—2irjlN\  ,  j  -  V— I.  The  sequence,  uin),  can  be  found  from  the  inverse  transfor¬ 
mation 

uin)  =  £  u{k)  wsk"-  «  -  0.  1,  2 . N  -  1.  (A5) 

Let  us  define  to  be  an  W  x  Af  matrix  such  that 

fiv  -  </*„>;  /™  =  m,n  =1,2 . /V.  (A6) 

The  matrix  F*  wiH  be  called  the  Mh  order  DFS  matrix  because  we  can  rewrite  Eq.  (A4)  in  matrix  nota¬ 
tion  as 

U  -  Fnu  (A7) 

where  U  -  (UiO),  U(  1) .  UiN  -  l))r,  u  -  (w(0),  «(1) . wUV-  l))r,  and  F  denotes  tran¬ 

spose.  The  DFS  matrix  has  the  property  that 

n-  (A8) 

This  property  can  be  shown  by  rewriting  the  inverse  DFS  transformation,  Eq.  (A5),  in  matrix  notation 
and  comparing  this  to  Eq.  (A7). 
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Let  us  define  a  periodic  sequence  y(/t),  n  -  0,  1,  2 . N  -  1  such  that 

on  n  -  0,  l . M  -  1 

y(n)  “  n  =  M,  M  +  \ . AT  -  1. 

It  can  be  shown  that  K(/c),  k  -  0,  1 . N  -  1  are  real. 


(A9) 


We  can  now  rewrite  Eq.  (A2)  in  matrix  notation  and  show  that  z  =  Ax  or  equivalently 

x  =  A~l  z  (A10) 


where  A  is  an  N  x  N  UDHTM  defined  by  Eq.  (Al),  z  -  (z(0),  z(N  —  l))r,  and  x  =  (x(0),  . . .  , 
x(N  —  l))r.  We  can  also  write  Eq.  (A3)  in  matrix  notation  as 

Z=A*  (All) 


where  A  is  a  diagonal  matrix  with  real  diagonal  elements  Kkk  -  Y(k  -  1),  k  =  . . .  ,  A,  Z  -  (Z(0), 

. . .  ,  Z(N  —  l))r,  and  X  =  (A'(O),  . ..  ,  X(N  —  l))r.  However,  we  know  that  Z  =  FAz  and  X  =  Fsx , 
so  that  Eq.  (All)  can  be  rewritten  as 

Fsz  =  \Fsx.  (A12) 

If  we  solve  for  x  in  Eq.  (A  12)  and  use  Eq.  (A8),  we  find  that 

X  =  ±:F;.\-'Fsz.  <  A 1 3 ) 


Subtracting  Eq.  (A13)  from  Eq.  (A10),  we  see  that 


0  = 


-i 


(AI4) 


Since  z  can  be  chosen  arbitrarily,  this  implies  that 

A-'-jj/Hi-’F*. 


(A15) 


We  summarize  our  result  by  the  following  theorem: 

Theorem-.  If  A  is  an  N  x  N  UDHTM,  then  A  can  be  written  in  the  form 

where  FN  is  the  Mh  order  DSF  matrix  and  A  is  an  N  x  N  diagonal  matrix  with  real  elements.  In  addi¬ 
tion,  if  A  is  nonsingular,  then  A~x  has  the  form 

F'"  A_l  F»- 


We  now  prove  the  following  theorem: 

Theorem:  If  A  is  a  nonsingular  UDHTM,  then  A~]  is  a  UDHTM. 
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Proof: 


Let  us  derive  an  individual  element  of  the  matrix  A  by  using  Eq.  (A1S).  By  direct  calculation, 
it  can  be  shown  that  if  A~x  -  ( amn )  m.n  -  1.  2 . N,  then 

^Z^t  WkH-m)k.  (A16) 


*-o 


Thus 


We  show  that  A'1  is  hermitian  by  using  the  fact  that  Kk  k  =  1.  ...N  is  real  and  W7 1  - 


77  1  w}r-"k \ 

A- 0 


n  k- 0  ' 


i  -  I 

_  _L  V  \  - 1  u,-(n-mU 

N  L  A*  +  > 


*-0 


Also,  it  is  readily  shown  from  Eq.  (A16)  that  the  diagonal  elements  ( m~n >  are  real. 

We  use  the  form  of  Eq.  (A  16)  to  show  that  A  1  is  Toepiltz.  We  see  that  it  is  possible  to  write 

in  the  form  «m„  =  fim  „  for  all  m  and  «,  which  is  exactly  the  form  of  a  Toepiltz  matrix. 

We  show  that  A  ' 1  has  the  up-down  property  by  demonstrating  that  for  the  elements  in  the  first 

row  that 

aln  “  «l<\-»  +  2>’  (A17) 


We  do  this  as  follows: 


a  KA-n  +  2)  ~ 


77  W+" 

”  A- 0 


-  77  t  ,  (»rr-',+l)f 

ly  k-0 
|  JV-I 

_  _L  V  x  -1  u/tn-l)* 

k,  L  A*  +  l 

i.  _  n 


k-0 
“  a  Iff' 

Hence,  the  theorem  is  proved. 

We  see  from  Eq.  (A16)  that  in  order  to  find  the  elements  of  the  first  row  of  A~\  we  can  write 

“u"  77  (A18> 

”  (k-0 


However,  we  notice  that  the  form  of  Eq.  (A  18)  is  that  of  a  DFS  (see  Eq.  (A4))  except  for  a  scalar  fac¬ 
tor  of  1/A.  Hence,  to  generate  the  first  row  of  A~\  we  merely  find  the  Mh  order  DFS  of  the  sequence 
X f *,  XjT ....  X  v1  and  divide  all  elements  of  the  DFS  by  N.  Therefore,  since  the  first  row  of  a  hermi¬ 
tian  Toepiltz  matrix  specifies  the  entire  matrix,  we  have  found  a  simple  method  of  generating  the 
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inverse  of  A~x.  Firstly,  we  generate  X2 . X*  by  calculating  the  Mh  order  DFS  of  the  sequence 

. •  ■  ■.«!*•  Secondly,  the  first  row  of  A~x  is  found  by  calculating  the  Mh  order  DFS 

of  the  sequence  Xf1,  X^1 . X*1  and  then  dividing  all  elements  of  the  DFS  by  N.  Finally,  because 

A~l  is  a  hermitian  Toepiltz  matrix,  all  other  elements  of  the  matrix  are  specified  by  elements  of  the 
first  row. 
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Appendix  B 

PROGRAM  LISTING  OF  THE  MATRIX  INVERSION  ALGORITHM 


C 

c 


c 

c 

c 

c 


SUBROUTINE  TOSPLZ  ( H,Gft  EGA ,  .1 ,  N) 


TFIS  SUDItOUTINF  FINDS  Tn?  xuVEaSZ  Ok  A  jxV-n  ..  jJSINGULAF 
KEPMITIAN  TOSPILT  Z  DAT  FIX  «.a„&£ 

H  =  IriS  HEFalTIAN  TCTrlETr.  DAI  r.lX 
(NOTE  THE  3  _  AGC  N  A  L  EL-ft-NIS  Mool 


C 


C 

c 


0  M  E  G .Y=  Iti  E  MATBIX  INVERSE  of  u 

K*IH3  PC*  OF  COLUMN  DxaEM  ION  Or  ;i 
N=2M— 1 


C 

r 


C 


c 


c 


liiE  ALGOPITHh  ;iAi  FAIL  IF  Am  I.«-  -m _DI .il -  a.t.A  I.iAT  IS 
JS  r  D  IN  CALC Ox Ax  _  NO  Till.  »  A  V  „  13  £  x£  S  _i«  _>  J  _/i_  .  ..r'  T  H  »  3 

OCUFRS  I  HE  MESSAGE  "ALGOFIIi.I  r'nlLs"  lx 


*  *  > 


IMPE  IC  IT  CO  MI  II  X“d  (x>— H  t  0 ,  *— x. j 

D  IM  ENSIGN  rf  £!,.•»)  ,  T  (N  ,A)  ,r  (H,  NJ  ,  IP  i  j  ,  1 1  ,  j14  JG  A  (  m 
1  w  v N  f  N )  /  v  1 1  v  ft  »  ft)  »  v  ^  2  (ft  ,  ri  J  i  i/x  i  v M  /  I )  ,  a  u  j  k  j  ,  a  )  ,  ( N 

1  THETA  (N  fS) 

t  AT  A  PI/ 3.  1413^29/,  A.  EE  r.  Jr/.  OoOoal  / 

IN  III ALIZ £  MATFIX  CONSTANTS 
INITIALIZE  I  MAT  FIX 
DO  300  K=1,K 
T(K,M)=H(K,1) 


300  CONTINUE 

M 1  =  M ♦ 1 

DO  400  K=  X1  f  N 
T(K,ft)=CONJG(ti  (S*2-K,1j) 
400  CONTINUE 


ft 1  =  ft—  1 

FINE  SUCCESSIVE  DPS 
DO  500 
M  M  =  M  ♦  1  —ft  ft  ft 
NN=2.*MM-1 


COMPUTE  DFS  MAT? IX  OF  0FD3E  NN 
Al=-2. *?I/NK 
DO  oOO  K= 1 , NN 
DO  700  L=1,NN’ 

AC=COS ( (K-1) *  (L-1)  *A1) 

AS=S IN  (  (L— 1 )  *  A1 ) 

F(K,L)  -CM PL X  (AC  , A 5) 
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7  ju  c  i  Nil  SUE 
oUO  CONTINUE 

JO  3000  h-= 1  ,  N  N 
3000  TF  (K,  1)  =T  (K 
C  FIND  Jr'S  OF  TP 

DO  32UJ  K=i,NN 
T  C  K  ,  1 )  =CMFLX(u.  ,0.) 

DO  33uO  1  =  1,  NN 

-«3JU  ?  (f>,1  )=F  (*,’)  ♦?  (K,  I)  *TF  U,  1j 
-lOu  CONTINUE 

JO  8oU  K=  1,  NS 

r =r .EM.:.  :k,^) 

-  U  )=■!./  (t-S*?.  IK,  -  ) ) 

30  0  CO  ft  I  IN  U  E 

DC  j40*j  *»=  1,  NS 

(K,‘ )  = c  k  ?  l  a  (0.  ,0.  j 
JO  350  0  1=1,  :u. 

3  5  j  u  :?  (K,1)  =  IF  ;  ft,1!  *F  (K,L)  *s(w  1} 
A  =  A oS  (?r.\L  (I  P  (h,-)  )  ) 

IF  JA.LT. AF'r.t.O?)  "?YP.*:  **00 
“•JO  F  3r  n  AT  ("  A  ,  '  *  LGCBI?  FAILS' 
if  (A.LT.At'I-  r  o:-)  SETUP., 

?<*J0  CONTINUE 

rpp  ,D*.3rALjTPr  D) 

-  r  »  FM)  =  TP  tl  ,1j 
DO  Sou  f>=i,NN 

:?jk,‘')=ccaoo(tp:k,1)) 

HKf:-.M)  =  Tr(h,’) 

300  CONTINUE 
«K1=»V.«-1 
00  1000  K=  1,.V.E  1 
I(K,:'iK‘»)=T?CA,‘') 

1000  CONTINUE 

ns  1=2.  *aa  i-l 

DO  HOO  K  =  M ,  N 1 
T  (K  ,;•]«'»  )=CONJG  (TP  (NN1»2-K,  1)  J 
HOO  CONTINUE 
500  CONIINUE 

C  CO fj ? U I  E  INVEFSE  MATEIX 

OilSGAJI  ,1)=1./T  (i,2> 

LO  1200  Mt!=2,a 
NN  =  2.  *AA- 1 
C  FIND  J  KA"FIX 

o(i  ,‘t>  =  r  (1  ,ea> 

DC  13o0  K  =  2,  NN 

1  300  U  T  ,K)  =CONJG  JT  JK,3fi)  ) 

CO  IwOO  N=2,.NN 
'CwO  w  JK,1)  =CCNJG(C  J1,K)  ) 

LO  1500  J=2 ,N  N 
DC  1600  1=2, NN 
W  (I  ,J)  =tt(I-%J-1) 

IcJo  CONTINUE 
*■500  CONTINUE 
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C  FIND  Oil  ,^2 
aal=A8-i 
MMPl=feM*1 
DO  3600  K  =  1 ,  M  K 
DO  3700  L=1  ,rttt 
3700  g11(K,L)*i2(K«L) 

3600  CONTINUE 

DO  3800  K=i  ,fiM 
DO  3900  L=1,MKl 
3900  g?2  (K,L)=0(K,*M*L) 

3800  CONTINUE 

DO  4000  K=  1,  MM  1 
DO  4100  L=i ,MM 
4100  w21  (K,L)  =U  (KM+K.L) 

40^0  CONTINUE 
C  COMPUTE  INVERSE 

DO  4200  1=1, *M1 

DO  4300  0=" ,MM 

H  80  (I, J) =  CMPL  X  (0. ,0.) 

DO  4400  K  =  1,MM1 

44  OvJ  ?HC  (I  ,  J)  =  FUG  (I  ,  J)  ♦CKEUA  JI,  iv)  *k21  [K,  0  > 
4300  CONTINUE 
4200  CONTINUE 

DO  45 00  1=1  ,  MM 
DO  4600  J  =  1<M!' 

THETA  (I,J)=COPLX  [0.  ,0.) 

CO  4700  K=1 ,KMl 

4  7uO  T  HIT?  ( I,  J)  =THiTA  JI ,  J)  ♦  ^1 2  il ,  K)  *FiiC((v,4| 
4600  CONTINUE 
4500  CONTINUE 

DO  4800  1=1, KM 
DO  4900  J=i  , MM 

4  90  0  OMEGA  (I,J)=g11  (I, Jj-TKETA  (I,J) 

4800  CONTINUE 
1200  CONTINUE 
HETUP.N 
END 


IS 


